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STABILITY OF PLANE POISEUILLE FLOW WITH HEAT TRANSFER 


by Merle C. Potter* and Edwin J. Graber, Jr. 

Lewis Research Center 

SUMMARY 

The influence of heat transfer on the stability of a plane Poiseuille flow in a channel 
is analytically investigated. The primary liquid flow is affected by the heat transfer 
through the variation in viscosity with temperature. Additional terms resulting from 
the viscosity gradient are included in the development of a modified Orr-Sommerfeld 
equation. It is the presence of these terms which leads to a prediction of more unstable 
flows, for without their inclusion, the flow is stabilized. The results show that a tem- 
perature difference between the walls is always destabilizing and in particular a temper- 
ature difference between the walls of 100° F (55. 6 K) leads to a reduction in the critical 
Reynolds number from 7800 to 4600. 


INTRODUCTION 

The linear stability characteristics of plane Poiseuille flow of a liquid between 
plates of fixed but different temperatures are investigated by the method of small dis- 
turbances. Particularly interesting is the influence of a temperature gradient on the 
minimum point (critical point) of the neutral stability curve. Since the viscosity of li- 
quids is temperature dependent, the imposed temperature gradient creates a viscosity 
variation in the flowing fluid. This gradient subsequently causes additional terms to ap- 
pear in the governing equation thereby altering the stability characteristics of the flow. 
To establish the importance of these additional terms, the stability equation is first 
solved with their inclusion and then with their omission. The primary velocity distribu- 
tion used in both cases is the skew symmetric profile which results from the inclusion 
of the viscosity variation in the solution of the primary flow momentum equation. 

Thomas (ref. 1), along with many others (e. g. , ref. 2), has analytically deter- 
mined the critical Reynolds number for isothermal Poiseuille flow in a channel. Thus 

♦Professor of Mechanical Engineering, Michigan State University, East Lansing, 
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in the present study the limiting case of equal wall temperatures (AT = 0) can be com- 
pared with Thomas' results. 

Potter and Smith (ref. 3) introduced a skew- symmetric primary velocity profile into 
the Orr-Sommerfeld equation (stability equation neglecting a viscosity variation) for iso- 
thermal Poiseuille flow and established that a double critical point exists for this case. 
The profile used in their analysis, however, was taken from experimental results and 
was not analytically shown to correspond to the stability equation used. It was also con- 
cluded by Potter and Smith that a skew- symmetric primary velocity profile creates a 
more stable flow situation than does the symmetric velocity distribution (that is, the 
Reynolds number increases as the flow deviates from the symmetric velocity profile). 

The present study was arranged to analytically produce a skew-symmetric velocity 
profile in water and then establish whether or not the corresponding stability equation 
yields a double critical point and results in a stabilizing effect as in Potter and Smith. 

The governing equation and corresponding boundary conditions form an eigenvalue 
problem, and the eigenvalues being the Reynolds number, the wave number, and the 
wave speed. It is the eigenvalues forming the neutral stability curves that are of pri- 
mary interest in this study. In particular, it investigates the influence of heating or 
cooling one wall on the minimum point of the neutral stability curve. The range of wall 
temperature differences considered is from AT = 0° F (0 K) to AT = 200° F (111. 1 K) 
(the symbols are defined in the appendix). 


ANALYSIS 
Primary Flow 

For steady two-dimensional flow of a viscous liquid flowing between heated parallel 
plates (see fig. 1), the x-component of the momentum equation reduces to 

(i) 

dY V dY/ dX 

Neglecting viscous dissipation (see ref. 4), the energy equation yields the linear temper- 
ature profile 


T - (T 2 - Tj)I + Tj (2) 

where Tg and Tj are the upper and lower plate temperatures, respectively, and d is 
the distance between the plates. 
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The viscosity relation used in this study is the conventional one for liquids, namely, 


Cj/T 

M — c 2^0 e 


(3) 


where c 2 and c^ are constants determined from viscosity-temperature curves and 
juq is the viscosity at the cool wall. 

The equations are nondimensionalized by using 


x 




u = 




~ a 
M = — 


Mq 


R = 


U n/> d 


h 0 



where U m is the average velocity. Equation (1) is, in dimensionless form, 


c 0 — <exp 
1 dy 


l (T 2 - T^y + Tj 


du l = R^P 

dy dx 


(4) 


The solution of this equation, which satisfies the no- slip boundary conditions at the 
plates, is 


u(y) 


R — 
dx 

2a 2 Cc 


(ay + b)(ay - b + A - l)exp (- — - — + b(b + 1 - A) exp (- - 


\ ay + b 


V b 


+ (A - 1 - 2b) 


El - 


1 


b / \ay + b. 


(5) 


where 


a = 



T 


1 


C 1 
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b =■ 


A = 


(a + b)(b + 1 - a)exp(- 

V a + b. 


- b(b + l)exp 


(- 5 ) 


+ (2b + 1) 


Ef 

.b / Va + b/J 


(a + b)exp( — 

\ a + b 


b exp (- -'j + E^— ^ - E 

b / Vb/ \a + b 


and E(x) = J' e"^(dt/t). In equation (5) R(dp/dx) can be eliminated by the continuity 

X 

/' 1 9 

relation u(y)dy = 1. As T 2 — TV R(dp/dx) — -12 and u — -6(y - y). Several 
of the velocity profiles are plotted in figure 1. The values used for the constants c^ 
and Cg were 3233° R (1796. 1 K) and 1/445.9, respectively. 


Stability Equation 


Let the stream function for an infinitesimal disturbance (a Tollmein-Schlichting 
wave) be given by 


+ = <P( y)e ia < x - cl > 


( 6 ) 


where a is the wave number and c = c . + ic-, the complex wave speed. It can be 
shown that Squire's Theorem is applicable; hence the use of only two-dimensional dis- 
turbances. Using u = u(y) + (d\p/dy) and v = -dip/dx as the perturbed velocity field in 
the Navier-Stokes equations results in, after nonlinear terms in dify/dx and d\p/dy are 
neglected and the pressure is eliminated by cross-differentiation, the modified Orr- 
Sommerfeld equation (ref. 5) 


(u - c ){(p" - a cp) - u "cp 


1±L {cp 

aR 


iv 0 2 ,, , 4 x 

la cp + a cp) 


1 

aR 


2 {cp'" 

dy 


a cp') + — — H {cp" + a‘ 

dy 2 


(v) 


in which R = U m pd//iQ. This equation differs from the usual Orr-Sommerfeld equations 
in that it contains the viscosity derivative terms. These terms are very significant for 
the study of liquids and must be retained. 


4 


The boundary conditions representing the no- slip conditions at the walls are 


cp(0) = cp'(O) = <p{ 1) = <p'(l) = 0 (8) 

The eigenvalue problem thus formed requires that 

c i (c r , a, R) = 0 (9) 

if no growth or decay of the disturbance is allowed. The minimum Reynolds number 
from the neutral- stability curve representing the previous relation is the critical Rey- 
nolds number sought in this study. It will be a function of the AT chosen between the 
two plates. The eigenfunction is not symmetric because of the skew- symmetric tem- 
perature imposed, hence one must use the whole channel in the analysis. Equation (7) 
with boundary conditions (8), the eigenvalue problem, must now be solved to yield the 
desired eigenvalues. A numerical technique, outlined in the following section, was 
chosen to solve the problem. 


Numerical Integration of the Modified Orr-Sommerfeld Equation 

The task of solving the modified Orr-Sommerfeld equation is made difficult by the 

_4 

presence of a very small coefficient of the highest order derivative of the order 10 
The asymptotic method developed by Heisenberg (ref. 6) and improved by Lin (refs. 

7 to 9) has been very popular but has certain limitations in addition to the difficulty in- 
volved by introducing the "critical" point. Numerical methods used to solve the Orr- 
Sommerfeld equation are proving to be quite successful. These methods generally re- 
quire multiple precision as did that used by Thomas (ref. 1). However, an initial value 
scheme, requiring only single precision, has been devised by Kaplan (ref. 10). This 
scheme, used by Reynolds and Potter (ref. 11), will now be outlined. 

Two linear independent solutions satisfying the lower wall boundary conditions must 
be combined to produce a solution satisfying the upper wall boundary conditions. One of 
these solutions grows very rapidly away from the wall, making it difficult to maintain 
two independent solutions during the numerical integration. To insure independency, 
Kaplan suggested suppressing the growing solution at each step in the calculations, that 
is, at each step a small multiple of the growing solution is subtracted from the behaved 
solution. The effect of the suppressed parts are later accounted for in the solutions. 
The two solutions thus formed are combined to satisfy the upper wall boundary condi- 
tions. 
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The algorithm employed in the numerical integration is part of a predictor- 
corrector algorithm for fourth-order equations. However, the linearity of the equations 
eliminates the need for a predictor, and hence only the corrector equations are re- 
quired. They are obtained by passing a third order polynomial through cp iv at four 
points, expressing the coefficients in terms of the values of the fourth derivative at the 
three backward points and the single forward point. This is then integrated to give cp ™ 
at the forward point in terms of the unknown fourth derivative, again to get <p" , and so 
forth. The resulting expressions for cp and its derivatives are 

2 3 4 

<t x = „„ + y 0 a + v'o j- + <P’ 0 " f + ~ W + 214 <f o V - 32 ^- V i + «*>*) 

"•] - ■’<] * - ! o * 1 "o' + f 2(l (’ i, ”'r + nn ‘ f » - 2 i « > - 1 i + 4< p'X) 

2 

tt 1 1 ttt \ A / 00 iv lr71 iv oc iv a „ iv \ 
cp ^ — cp q + cp q A + i 3Scp ^ + 171c/? q - 3 6 ^ + 7 cp 

360 ' 

<p'l = <Pg' +— (vi V + 19 ^q V - + <^ 2 ) 

2 4 

where A is the mesh size and cp = <p(y n ). These expressions are substituted into the 
modified Orr-Sommerfeld equation yielding an equation in which cp™ is the only un- 
known. A similar scheme, based on a two-point fit, is used to start the calculation. 

To use the previously outlined method, the eigenvalues must be known. If the eigen- 
values are known only approximately, the upper wall boundary conditions cannot be 
satisfied exactly. To find the exact eigenvalues, the following iterative scheme is used. 

Let (p T,T (l) = 0 and choose cp{l) as a test function E. The object is to make E zero 
_ <1 

(10 will suffice). Fixing the wave speed c(c^ = 0), choose Aa and AR such that the 
positive-definite quantity EE is minimized (E is the complex conjugate of E). 



where ?E/3o and 3E/PR are found by using 0. 01 percent changes in a and R. Con- 
vergence to the exact values is quite rapid requiring only three or four passes. If a 
Control Data 3600 computer were used where each pass required 8 seconds, a neutral 
stability curve could be found in approximately 10 minutes if initial guesses were made 
carefully. 
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RESULTS AND CONCLUSIONS 


The primary velocity profiles corresponding to an imposed wall temperature differ- 
ence of from AT = 0° F (0 K) to AT = 200° F (111. 1 K) are given in dimensionless 
form in figure 1. Since the viscosity of liquids decreases with increasing temperature, 
the resultant velocity profiles became skewed as | AT | increases; the maximum veloc- 
ity is shifted towards the hotter wall. 

Figure 2 shows that, if the viscosity gradient terms are neglected in the stability 
analysis, the resultant neutral stability curve is shifted to the right as | AT | increases, 
thus indicating a more stable flow situation. This agrees with the results of Potter and 
Smith (ref. 3). In contrast, when the terms are accounted for, the flow becomes de- 
stabilized as | AT | increases. It should also be noted that a double critical point did 
not appear in the neutral stability curve when the viscosity gradient terms were included 
in the equation. However, neglecting the viscous terms led to the appearance of an in- 
flection point (or possibly a double critical point for other AT's) at a Reynolds number 
of approximately 32 000 for AT = 100° F (55. 6 K) (fig. 2). Hence, even though the vis- 
cosity gradient terms are small, their inclusion is extremely important. 

Figure 3 plots the variation of the critical Reynolds number with wall temperature 
difference AT. It can be seen that the temperature difference between the plates has a 
definite destabilizing effect; there is a 50 percent reduction in critical Reynolds number 
resulting from a 140° F (77. 8 K) temperature difference. 

The curves for AT = 0° F (0 K) (figs. 2 and 4) represent isothermal Poiseuille flow 
and agree with the accepted results of Thomas (ref. 1). Eigenfunctions are also plotted 
in figure 5 for AT = 0° F (0 K) and for AT - 200° F (111. 1 K). 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, July 3, 1970, 

720-03. 
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APPENDIX - SYMBOLS 


A 

constant 

U m 

average axial liquid velocity 

a 

constant, (Tg - T^/c^ 

u 

dimensionless axial velocity, U/Uj 

b 

constant, Tj/c^ 

u 

perturbed axial velocity 

c 

complex wave speed, c r + ic^ 

V 

perturbed transverse velocity 

c l> c 2 

constants 

X 

axial coordinate 

d 

distance between plates 

X 

dimensionless axial coordinate, 

P 

static pressure of liquid 


X/d 



Y 

transverse coordinate 

P 

dimensionless static pressure, 




P/puf n 

y 

dimensionless transverse coordi- 




nate, Y/d 

R 

Reynolds number, U pd/ p n 




m u 

a 

wave number 

R nr 

critical Reynolds number 



cr 


A 

mesh size 

T 

liquid temperature 





M 

viscosity 

AT 

wall temperature difference, 




T 2- T l 

M 

dimensionless viscosity, p/p q 

T 1 

lower wall temperature 

M 0 

viscosity at 

T 2 

upper wall temperature 

P 

density 

t 

time 


eigenfunction 



iL/ 

stream function 

U 

axial velocity of liquid 

T 
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Figure 1. - Primary velocity distribution. 
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Figure 2. - Neutral stability curves for various temperature differences AT. 



Figure 3. - Variation of critical Reynolds number with temperature difference AT. 
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Figure 4. - Variation of wave speed with Reynolds number for various temperature differ 
ences AT. 
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Figure 5. - Eigenfunctions at critical point 
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